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The description of a stellar system as a continuous fluid represents a convenient 
first approximation to stellar dynamics, and its derivation from the kinetic theory 
is standard. The challenge lies in providing adequate closure approximations for 
the higher-order moments of the phase-space density function that appear in the 
fluid dynamical equations. Such closure approximations may be found using rep- 
resentations of the phase-space density as embodied in the kinetic theory. In the 
classic approach of Chapman and Enskog, one is led to the Navier-Stokes equa- 
tions, which are known to be inaccurate when the mean free paths of particles are 
long, as they are in many stellar systems. To improve on the fluid description, we 
derive here a modified closure relation using a Fokker-Planck collision operator. To 
illustrate the nature of our approximation, we apply it to the study of gravitational 
instability. The instability proceeds in a qualitative manner as given by the Navier- 
Stokes equations but, in our description, the damped modes are considerably closer 
to marginality, especially at small scales. 



I. A KINETIC EQUATION 

If we have a system of N stars, with N very large, and wish to study its large-scale 
dynamics, we have to choose the level of detail we can profitably treat. Even if we could 
know the positions and velocities of all N stars for all times, we would be mainly interested 
in the global properties that are implied by this information. For just such reasons, many 
investigators prefer to find an approach that leads directly to a macroscopic description of 
the dynamics. As Ogorodnikov (1965) has put it, "In order to exhibit more clearly the 
kinematics of highly rarefied media, and of stellar systems in particular, it is useful to make 
a comparison with the motion of a fluid." However, traditional methods for deriving fluid 
equations are effective only for media in which the mean free paths of constituent particles 
are very short compared to all macroscopic scales of interest. This condition is not met in 
many stellar systems and plasmas and so we here describe an approach that is effective for 
deriving fluid equations for rarefied media such as stellar systems. 

The first problem we must face is to decide what the kinetic description of a stellar dy- 
namical system ought to be. Since an A^-body description is not what we want to work with, 
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even if we could, since that approach would have us computing the complicated trajectory 
of the stellar system through a phase space with large dimension. So we go straight to 
the description of the system in the six-dimensional phase space whose coordinates are the 
spatial coordinates (x^) and the velocities (v*) of the N stars, where i = 1, 2, 3. 

A plot of the locations of each of the N stars in the six-dimensional phase space at some 
given time, would reveal a swarm of points whose detailed description would also be too 
complicated for us, at least in a first look at the problem. So instead, we concentrate on an 
ensemble mean of such a description and seek an equation for the density distribution of this 
mean. That equation, on which we base this work, is an evolution or continuity equation 
for the probability density in the six-dimensional phase space. 

Since we treat the stars as points, the true density F in the six-dimensional phase space 
is a summation of delta functions at suitable locations. This density is advected by a six- 
dimensional phase velocity that Hamilton's equations tell us is solenoidal. Hence the total 
time derivative of the density is 

DF = (1) 
where the comoving derivative in phase space is defined as 

■D := 9t + V • Vx + a ■ Vv (2) 

and the subscripts x and v on the gradient symbols indicate that they are gradients with 
respect to position and velocity respectively. As usual, v = x where the dot means total 
time derivative and the quantity a is the gravitational acceleration per unit mass; we assume 
that all the stars in the system have the same mass, m. The gravitational acceleration is 
given by the gradient of the gravitational potential per unit mass, which is a solution of 
Poisson's equation, 

A$ = 47rG' / Fd\. (3) 



Let /(x, V, t) be the ensemble mean of F. Then the total density will he F = f + f where 
/ represents the fluctuations about the ensemble mean; / will involve the same summation of 
delta functions as does F plus a smooth background distribution with negative mass density 
arranged so that the ensemble average of / is zero. We similarly split the gravitational 
potential $ into an ensemble mean part (p plus a fluctuating part (p. Then, if we take the 
ensemble average of (1), we obtain 

dtf + v Vx/ + (Vx0) • Vv/ = -ii^J) ■ Vv/) , (4) 

where 

A(f) = 47rG /" / dh. (5) 



The terms on the left hand side of (4) describe how the mean phase density, or distribution 
function, /, streams through the single-particle phase space. The right side represents the 
mean influence on the evolution of / exerted by the average of the fluctuation-interaction 
term. The latter represents the self-interactions of the system caused by fluctuating effects 
and may be thought of as representing the influence of collective modes (such as waves or 
quasiparticlcs) on the motions of the individual particles. 

It is typical that the chance of close approaches of two stars in many stellar systems is 
small. Because of this, one commonly made approximation is to neglect the right side of (4) 
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completely and so work with what is called the coUisionless Boltzmann (or Vlasov) equation. 
More reasonably perhaps, one may try to give an expression for the way the self-interaction 
term affects the flow of the phase density through the phase space. 

As in Boltzmann theory, we shall suppose that the right side of (4) may be expressed as 
a functional of / itself so that the kinetic equation is deterministic. That is, we assume that 
the kinetic equation takes the form 

dtf + y V./ + (Vx0) • Vv/ = e[/] (6) 

where C[-] may be called a collision term in keeping with the terminology of kinetic theory 
even though it does not arise from direct binary collisions. This approach may be acceptable 
because it appears that the fluid description that we seek is not very sensitive to the details 
of the right side of (4). On the other hand, we must admit that this hope is founded on a 
very limited range of trial forms since the job of deriving the consequences of each form is 
laborious. Our aim here is to adopt one standard form for the collision term and use it to 
go on to a coarser description of the stellar system like that of fluid dynamics. 

The parallel to the Boltzmann theory has been used to good effect in the study of plasma 
physics, as in the work of Rosenbluth et al. (1957). As Clemmow and Dougherty (1969) 
explain, those authors obtained their results by "expanding the Boltzmann collision oper- 
ator under the approximation that all the deflections are small angle and cutting off the 
impact parameter at about the Debye length. ... At first sight the success of that method 
is surprising, as any treatment dealing with binary collisions would seem to be discredited. 
The physical reason for the agreement is that, for the majority of particles, there is little dif- 
ference between a succession of numerous small-angle collisions (regarded as instantaneous 
and occurring at random) and the stochastic deflections due to the presence of many nearby 
particles continually exerting weak forces. These two pictures of the dynamics are of course 
represented respectively by the Boltzmann and the Fokker-Planck" approaches. 

Similar thoughts have been expressed in the context of stellar dynamics, most recently 
by Griv et al. (2001) and by kinetic theorists generally. As E. G. D. Cohen reports (1997), 
"when Academician Bogolubov and I discussed the nature of kinetic equations, he mentioned 
a discussion he had had with Professor A. Vlasov, where they had agreed that: Yes, in 
first approximation the kinetic equations for gases with strong short-range forces (i.e. the 
Boltzmann equation) and for gases with long-range forces (i.e. the Vlasov equation) differ, 
but in higher approximations they will become more and more similar. How right they were." 
A formal theory to buttress these remarks would be very comforting, but though we do not 
have one we shall adopt the point of view that the Fokker-Plank terms capture the essence of 
interactions of the stars in the system. Having thus supported our approach by the appeal to 
authority, we turn to the main purpose of this work, the derivation of continuum mechanics 
from the microscopic theory in a way that is not severely restricted to the case of short mean 
free paths. 

II. THE COLLISION TERM 

In the spirit of standard kinetic theories, we shall suppose that the effect of the coUision 
term is to drive the system toward a local equilibrium, though the correct equilibrium of a 
stellar system is not known on purely theoretical grounds. The tendency to approach an 
equilibrium seems not even to require a collision term of the usual kind since the violent 
relaxation described by Lynden-Bell apparently can do the job. Nevertheless, we shall 
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proceed in terms of the kinetic theory under the assumption that the spreading of the phase 
fluid through the phase space may be effected by a collision term. Furthermore, though the 
gravitational force has long range in physical space, we shall presume that this spreading 
takes the form of a diffusion of / through velocity space, that is, by the agency of a Fokker- 
Planck form of the collision term. (This may not be completely unfounded since there seems 
to exist a form of gravitational shielding (Spiegel, 1998) that may support the idea of local 
behaviour.) In this spirit, we write 



(7) 



The coefficients and B'^^ are generally functions of (x, v, t) and they may also be function- 
als of /. In this discussion, the choice of these coefficients in the Fokker-Planck description 
is adapted to the equilibrium that is expected or assumed. This equilibrium is a local one 
that satisfies the condition C[/o] = 0. 

Our goal is to find equations that govern the dynamics of the macroscopic properties of 
the fluid embodied in the density, temperature and velocity. These are defined as: 

Mass density P •= / i^fd^v, (8) 



Mean velocity u := -J mvfd^v, (9) 



Temperature T := / c^f dh, (10) 

3Rp J 

where the peculiar velocity is 

c(x,v,i) := v-u(x,i), (11) 

and R — k/m, k being Boltzmann's constant. 

As to the nature of C[/], we shall design it so that it produces what may be the simplest 
plausible equilibrium, namely the Maxwell-Boltzmann distribution 

= mi2nRTr/2 (-^) " (^2) 

Since T, p and u generally depend on x and this is a local equilibrium and we choose 
(Clemmow and Dougherty, 1969) 

A' = -T-^{v' - u'), B'^ = 2t-^RT5'\ (13) 

so that C[/o] = 0. We assume that the mean- free-time r is a constant so that the Fokker- 
Planck operator is linear in /. 

The collision term adopted here ensures the conservation of mass, momentum and energy 
in the system. This is reflected in the property 

V^"e[/]d3^; = 0, a = 0,...,4, (14) 



where 



V'" = m(l,v, 1^2). 



(15) 
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Thus we neglect the possible effects of dissipative processes and of evaporation of stars from 
the system. 

The macroscopic quantities (8)-(10) enter the equilibrium state (12) about which we are 
expanding. To ensure that the same macroscopic quantities that follow from / are those 
that determine fo, we impose a consistency requirement known as the matching conditions, 

Jrfd\^J^P"fod\, a^O,...,A. (16) 



III. FLUID EQUATIONS 

When we multiply the kinetic equation (6) by the coUisional invariants (15) and integrate 
over V, the right-hand side does not contribute to the outcome, and we are left with 

dtp + V ■ipu) = 0, (17) 
9tu + u- Vu = -p-^V-P- V0, (18) 
^pR {dtT+u ■ VT) = -P : Vu - V • q, (19) 

where V means Vx. Here the pressure tensor P and heat flux q are defined as 

F:^Jmccfdh, c^:^ j \mc^cf d?v. (20) 

We see that the form of the macroscopic equations is just that of the usual fluid equations. 
This result is independent of the rarity of the medium. The usefulness of these equations 
depends entirely on how well we can prescribe the higher-order moments P and q. A standard 
way to proceed is to solve (6) approximately for /. Wc shall follow this route also, but will 
deviate from the normally used prescription at a certain point. 

We let / = /o r/i . . . and look flrst at order r°. We flnd that e[/o] = 0, and the 
solution /o is the Maxwell-Boltzmann equilibrium (12). Prom (20), we see that Pq =pl 
and qo = 0, where the scalar pressure is given by p :— pRT. If we stop at this order. 
Equations (17)-(19) arc then the Euler equations for an ideal fluid. 

At order r^, it is convenient to factor out the Maxwell-Boltzmann solution from /i and 
write the equation to be solved as 

£A = 2)ln/o, (21) 

where A := /i//o and 

Lf := -c ■ Ve/ + i?TVe ■ (Ve/). (22) 
The right-hand side of (21) may be written out as 



2 

c , n 
+ In- 



2RT (27ri?r)3/2 
We note that 



^ f=-^c+fTT^-|)^lnT + Xilnp. (23) 



RT \2RT 2 



Dc = - V0 - Du/Dt - c • Vu, where D / Dt := dt + u-V . 



(24) 
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The operator £ maps polynomials in c to polynomials of the same degree so, given the 
form of (23), we may seek a solution /i to (21) as a cubic in c. We write 



/i = a + bi c* + +Cij cV + Dijk dc'c^ 



(25) 



where a, b, c and are functions of x and t, and symmetric in their indices; repeated indices 
are summed. Inserting (25) into the left-hand side of (21), we obtain 



(26) 



There is no a term because it is annihilated by the Fokker-Planck collision operator. We 
now equate coefficients of c between (26) and (23), and find 



1 



i8i?r 



{5ij V^k In T + 5ik V^i In T + 5kj V^i In T) , 



iRT 



Dt 



1 1 

bi = — V^i InT 

6 RT 



Dt 



- V^i In p. 



(27) 
(28) 
(29) 



Since a is still unspecified, we may use it to satisfy the matching condition (16) J /i d^v = 0, 
resulting from mass conservation. Only the terms even in c contribute, and we find 



/ 



hd?v^n{RTcu + o)^^, 



(30) 



which allows us to solve for a in terms of the trace of c. 

The pressure tensor and heat fiux (20) are then obtained from /i by performing straight- 
forward Gaussian integrals, and we get 



Pij = P Sij + 2pR Tcij , qk^lpRT{5bk + 21RT ^ee) ■ 



(31) 



The a term is absent from the pressure because we used the density matching condition (30). 
Prom (31) and (28), we find that the pressure tensor can be written 



P = pI-2//E-/x 



-^ + iV.u|I 



to first order in r, where the viscosity /i :— |pT, and 

^ij ■= \ ^xiUi Va^iMj - I V ■ U%) 

is the rate-of-strain tensor in traceless form. 

From (31), (27), and (29), to first order in r, the heat fiux is 



(32) 



(33) 



-77 VT- 377T 



V Inp -|- 



1 (Dw. 



- 



RT \Dt 



(34) 



where the thermal conductivity rj (5/6)p ri?. 
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These results differ from those of the usual Navier-Stokes equations for which ¥ = pi — 
2fiK and q = —rj'VT. To get some understanding of the import of the additional terms 
found here we introduce the specific entropy 

S^C,\n{pp-'/^), (35) 

where C„ := 3R/2 is the specific heat at constant volume. Since 



Dt~ 



D\nT . 



(36) 



we find that 



P = p(^l-^'5jl-2/.E + O(r^) , (37) 
q^-r)VT + 3{r]T/p)V ■T + 0{t''), (38) 

with T := (P — pi). If we put these results into S, we obtain 

on 

^ = -— ^ [T: Vu+ V-ql. (39) 
3p 

Hence S can be seen to be 0(t) and so the additional terms in the pressure tensor are O(r^). 
A similar argument may be made for the new terms in the heat flux. Though these terms 
do not appear in the conventional fiuid equations, they can be quite significant when the 
mean free paths are long. 

If we eliminate the entropy using 

p 3p ^ ^ 

we encounter the combination p{t) — ^Tp{t) which are the first two terms of a Taylor series 
of p in t/2. We may then write 

P = p(i- ir)I + f //V-uI-2//E + 0(t2). (41) 

We see that our procedure has taken account of the physical fact that the medium senses 
what particles were doing one collision time prior to the present time but is not yet aware 
of what they are doing at the present. 



IV. THE JEANS INSTABILITY 

As an application of the equations of motion (17)-(19), we will use them together with 
the pressure tensor and heat fiux derived in Section III to examine the Jeans criterion 
for gravitational instability. This instability, describing the gravitational collapse of a ho- 
mogeneous medium, was first investigated by Jeans (1929). He found that perturbations 
above a critical wavelength (the Jeans length) were unstable to gravitational collapse, but 
that shorter wavelengths were unaffected due to the large-scale nature of the gravitational 
force. The Jeans length is the ratio of the adiabatic sound speed as to the gravitational 
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frequency y/AnGp. Pacholczyk and Stodolkiewicz (1959) and Kato and Kumar (1960) in- 
vestigated the effect of viscosity and thermal conductivity on the instabihty and found that 
the collapse occurs above a critical wavelength given by the ratio of the isothermal sound 
speed qt to the gravitational frequency. That critical wavelength is slightly larger than 
in the ideal case, which involves the adiabatic sound speed. Their interpretation for the 
increase in critical wavelength is that temperature gradients are adverse to the collapse, and 
a nonzero thermal conductivity allows the smoothing out of these gradients through very 
slow displacements of the medium; the mode is thus isothermal. 

We will now investigate how the gravitational instability occurs in our set of equations. 
We take as our equilibrium a medium at rest and with uniform density po and tempera- 
ture To- We expand each fluid variable into an equilibrium piece and a small perturbation, 

p = Po(l + <^), T^To(l + e), pRT ^po(l + w), w^^p + e. (42) 

The gravitational potential 4> in (18) is obtained from the Poisson equation 

A0 = 477^(^-^0), (43) 

where we choose — po as a background "neutralising" density in order to have a proper 
uniform equilibrium about which to expand; the equilibrium velocity then vanishes. The 
density —po is a repulsion term and may be regarded as a Newtonian analogue of Einstein's 
cosmo logical constant. To leave it out as Jeans did and jump straight to the linearised 
equation (47) below is expedient but questionable. 

The hnearised equations of motion (17)-(19) and Poisson equation (43) are 

dtip + V ■u = 0, (44) 

podtu + V ■ P = -poV0, (45) 

^PodtO + Po V ■ u + V • q - 0, (46) 

A4> = iirGpo ifi. (47) 

The pressure tensor P is given by (32) and the heat flux q by (34). We take the divergence 
of the velocity equation (45) and use the continuity equation (44) to eliminate V • u, 

Pod^ip-VV -.F^AnGpl^, (48) 

where we also used the Poisson equation (47) to eliminate (j). We then need to take two 
divergences of the linearised pressure tensor, 

VV : P = Po A((p + e)+ p(29tA(p - dtAe) . (49) 

We introduce the isothermal sound speed qt, the kinematic viscosity u, and the thermal 
diffusivity k, through 

ar — — , — — , « — —TT = (^0) 

Po Po Po '^p 2^P0 

Then, on inserting (49) into (48), we obtain 

9> = 4 ^{V + ^) + i^(29tA(p - dt^e) + Ik] a\ (p, (51) 
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where the Jeans wavenumber is given by 

^2 _ 47rGpo ^ 3 AnGpo ^^^^ 

Next, we take the divergence of the hnearised heat flux, 

V • q = -77^0 - STyToa^;^ [5* V • u + A0 + a^A((^ + e)] , (53) 

and insert this into the temperature equation (46), 

dtd - ptcp - Ik [Aa'^Ae - 3 {d'^(p - ^kj a^cp- a'^Acp)] = 0, (54) 

where again we ehminated V • u using the continuity equation (44), and used the defini- 
tion (50) of the thermal diffusivity k. 

Equations (51) and (54) are the equations required to derive a dispersion relation for 
the gravitational instability. It is convenient to use the viscous time fi/po as unit of time 
and aTlJ>/po as unit of length. Recycling the same symbols for the dimensionless quantities 
turns equations (51) and (54) into 

{d^ - A - 2dtA - |A:2) ip+{dt- I) AO = 0, (55) 
{3dt -fA)e+ {-2dt + 10 {d^ - Ik^ - A)) ip = 0. (56) 

On letting ip,9 exp{ikx + jt), we find the dispersion relation 

(3 + 2e) ^^ + ^ej'+[fk' + {i-f fc2) k' - e,] j + k^-i e^) e = o, (57) 

which may be compared to the expression obtained from Navier-Stokes by Kato & Ku- 
mar (1960), 

h' + f,k^l'+{lk' + k'-k])^ + l{k'-lk])e^^, (58) 

with the Fokker-Planck values for the viscosity and thermal diffusivity inserted into their 
result. In each case the system is marginally stable with 7 = at A;^ = (5/3)/cj, and is 
damped for larger k. This is illustrated in Fig. 1. For /c = 0, both dispersion relations 
predict a growth rate of 7(/c = 0) = >/5 /3 /c j in dimensionless units, which with dimensions 
is ^/^nGpQ. This is consistent with the fact that dissipation is unimportant at large scales, 
so the growth rate at /c = involves only the gravitational time. 

The asymptotic growth rate as /c — > 00 is —3/4 for Navier-Stokes and —1/11 for our sys- 
tem, independent of k. (Multiply by ^/pq to recover dimensions.) Thus, at large wavenum- 
bers, the modes tend to be uniformly damped, both in our case and for Navier-Stokes. This 
is because at large k the fluid behaves like a Stokes flow, where we can ignore the inertial 
and gravitational terms completely, and the balance is between the Laplacian of the pres- 
sure and the viscosity, which have the same number of spatial derivatives; hence the lack of 
dependence on k. For our case (the equations of Section III) there is a contribution from 
the new terms that shifts the damping rate considerably closer to marginality. 

Coexisting with the real root associated with the instability, there is also a pair of un- 
conditionally damped roots. The real part of these roots is plotted in Fig. 2. At /c = one 
of these roots is marginal (real part of growth rate equal to zero) and is an isopycnal mode 
(constant density). However, this mode is never destabilised and its growth rate immediately 
decreases as k increases away from zero. 
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FIG. 1 Growth rate 7 as a function of wavenumber k, for Navier-Stokes (dashed hne) and the 
equations derived in Section III (sohd). The growth rate of this mode is real and, in both cases, 
damping sets in above the isothermal Jeans wavenumber (5/3) A;j (here we have taken kj = 1 in 
dimensionless units). 

For small k the two damped roots are distinct and real, but they come together at 
larger k and become a complex-conjugate pair with nonzero imaginary part (not plotted), 
indicating oscillatory behaviour. The Navier-Stokes case (dashed line) is seen to have heavily 
damped complex roots at large wavenumber.^ But for our system of equation (solid line), 
the growth rate actually increases a little, looking as though it may be headed for a Hopf 
bifurcation {i.e., overst ability, where the real part of the growth rate becomes positive at 
nonzero imaginary part), before leveling off at an asymptotic value of the damping rate 
given by —53/165 ~ —0.3212. Indeed, it can be shown that there is no Hopf bifurcation for 
any realisable parameter values in our equations, but the fact that the complex modes are 
somewhat "destabilised" by the new terms is intriguing (this is also true to a lesser extent 
for the real mode described above). This destabilisation has its source in the k'^ coefficient 
of the 7^ term in the dispersion relation (57), which is not present in Navier-Stokes, but its 
physical significance is not yet apparent to us. 



V. CONCLUSION 

The basic approach in this as in other derivations of fluid equations from kinetic theory 
is to write the general moment equations (17)-(19). These are the fluid equations and, to 
complete them, we need closure relations for the pressure tensor and heat flux. This is an 
issue astronomers are familiar with from the study of radiative transfer. For the purpose, we 
could invent a phenomenological approximation as Eddington did in radiative transfer or we 



In the Navier-Stokes system the complex roots come together again at very large k, but the growth rate 
continues to decrease with k. 
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FIG. 2 Real part of the two complex roots 7 as a function of wavenumber fc, for Navier-Stokes 
(dashed line) and the equations derived in Section III (solid). Large wavenumbers are much more 
strongly damped for Navier-Stokes than in the system derived here. 



may pursue approximate solutions of the kinetic theory as Hilbert did by expanding in the 
collision time. The Hilbert expansion was developed by Chapman and Enskog in deriving 
the Navier-Stokes equations (Uhlenbeck and Ford, 1963) and we have pursued that line as 
well following earlier work (Chen, 2000; Chen et ai, 2000, 2001) based on the relaxation 
model of kinetic theory (Bhatnager et ai, 1954; Wclandcr, 1954). However, in that latter 
work, as here, we depart from the Chapman-Enskog approach in an essential way in not 
using results from lower orders to to simplify the results in the current order. 

To express this idea in equations, let us consider what happens in general in such problems 
in the first order. Once we have expressed the one-particle distribution function as / = 
/o(l + T(fi) where r is the (small) collision time, we are led to an equation for ip in the form 

= K/o + 0(t) (59) 

where £ is the linearisation of the collision operator. In general, £ is self-adjoint, as it is in 
the present study. Then fit/j"' — implies J ■i/j'^JCifdv — and so we must have 

j V'"[2)/o]rfv = 0(r). (60) 

Since ip"' represents the collisionally invariant quantities, the fiuid equations to the current 
order may serve as the solvability condition (60). 

In the Chapman-Enskog procedure, this solvability condition is taken to be a lower order 
version of the fluid equations, here the Euler equation, and it is used to simplify the right 
side of (59). Then, the results are used in the general fluid equations. For both of these two 
conditions to be satisfied, we require r to be very small indeed. 

What we are doing here is to say that, to the first order, the fiuid equations themselves are 
a realisation of condition (60) and that it is redundant to apply the same condition twice. 
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once with 0(r) retained and once with it omitted, as one does in the Chapman-Enskog 
method. Rather, we simply use the full condition (60) as a compatibility condition. It is for 
this reason that, in the first order theory, we allow ourselves to differ from the Chapman- 
Enskog results by terms O(r^). In particular, we have for the trace of the pressure tensor 



This result differs from the exact trace (3p) by O(t^) and this, we suggest is allowed in a 

first order-theory. 

We may add that in comparing the results of this approach to experiments on ultrasound 
we find that they do better than the usual Navier-Stokes version. Here we have an interesting 
example of a dictum of J. B. Keller: "Two theories may have the same accuracy but different 
domains of validity." 

We regret that though, in honour of Douglas Cough's birthday, we have gone to second 
order in this approach (for only the relaxation model so far), we could not fit the derivations 
into the space we were allotted in this volume. So those results will have to be presented 
elsewhere. We are happy to report that, in that next order, the trace of our pressure tensor 
differs from the exact trace by 0(t'^). For now we must be content with mentioning that 
result and presenting our best wishes to Douglas on his birthday. 
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